Genome-wide meta-analysis of QTL for morphological related traits of flag leaf in bread wheat

Flag leaf is an important organ for photosynthesis of wheat plants, and a key factor affecting wheat yield. In this study, quantitative trait loci (QTL) for flag leaf morphological traits in wheat reported since 2010 were collected to investigate the genetic mechanism of these traits. Integration of 304 QTLs from various mapping populations into a high-density consensus map composed of various types of molecular markers as well as QTL meta-analysis discovered 55 meta-QTLs (MQTL) controlling morphological traits of flag leaves, of which 10 MQTLs were confirmed by GWAS. Four high-confidence MQTLs (MQTL-1, MQTL-11, MQTL-13, and MQTL-52) were screened out from 55 MQTLs, with an average confidence interval of 0.82 cM and a physical distance of 9.4 Mb, according to the definition of hcMQTL. Ten wheat orthologs from rice (7) and Arabidopsis (3) that regulated leaf angle, development and morphogenesis traits were identified in the hcMQTL region using comparative genomics, and were speculated to be potential candidate genes regulating flag leaf morphological traits in wheat. The results from this study provides valuable information for fine mapping and molecular markers assisted selection to improve morphological characters in wheat flag leaf.


Introduction
Wheat is one of the world's three major crops, providing approximately a quarter of food for human. The continuous increase of wheat yield is crucial to meet the challenge of increasing food consumption [1]. Increasing planting density by improving the plant architecture of wheat on limited land is an effective strategy to increase yield [2]. In crops, canopy leaves, especially flag leaf, are the main source of dry matter accumulation in the grain filling stage [3,4], and flag leaf provide 41-43% of carbohydrates for grain filling [5]. Therefore, optimizing the morphological structure of flag leaves is a suitable method to improve plant architecture, photosynthetic efficiency and yield. Wheat flag leaf morphological traits are quantitative traits influenced by many environmental factors and controlled by multiple genes [6][7][8]. Many genes and QTLs that control leaf size and angle have been reported in rice [9][10][11][12][13]. For example, the mutation of OsDWARF gene resulted in the defect of brassinosteroid synthesis, which led to the reduction of plant height and the upright leaves [14]. In maize, Ku et al. [15] detected a major QTL qLA2 controlling the angle of flag leaf on chromosome 2. Tian et al. [16] found that UPA1 and UPA2 genes can increase the planting density by regulating the leaf angle of plants, thus increasing the maize yield. QTLs for flag leaf length, width, area and angle have been identified on 21 chromosomes in wheat [17][18][19][20][21][22]. For example, Liu et al. [19] detected three major QTLs on 3D, 7B and 7D for flag leaf angle. Liu et al. [23] found that TaSPL8 regulated leaf development by influencing auxin signal and brassinolide biosynthesis pathway, and affected flag leaf angle in wheat. Wang et al. [24] introduced the chromosome 1P of the wild related species Agropyron cristatum into common wheat to significantly reduce plant height and leaf size, thereby improving plant architecture and achieving dense planting.
Currently, many QTLs for flag leaf morphological related traits in wheat have been identified in previous studies. In order to make more effective use of the QTL for flag leaf morphological traits in wheat breeding, and deeply understand the genetic mechanism underlying flag leaf morphological traits, it is necessary to comprehensively analyze these QTLs to identify stable major genetic loci in wheat. QTL meta-analysis has been shown to be an effective method for integrating QTLs from various experiments onto a consensus map, narrowing QTL confidence intervals, and identifying reliable and stable meta-QTLs (MQTL) [25]. This method has been widely used in different crops for various traits, such as nematode resistance in soybean [26], yield under drought conditions in rice [27], yield and quality traits in cotton [28], yield in maize [29], and yield, nitrogen use efficiency, quality traits, disease resistance and abiotic stress tolerance in wheat [30][31][32][33][34][35][36][37].
With the development of high-throughput SNP sequencing technology, QTL mapping for complex quantitative traits based on natural populations using genome-wide association studies ( [42,43]. These studies indicated that the QTL location information identified by GWAS can effectively verify important QTLs, so that key genomic regions and candidate genes controlling important quantitative traits can be mined.
To date, QTL meta-analysis for flag leaf morphological traits has not been reported in wheat. In this study, QTL meta-analysis was performed based on QTL for flag leaf morphological traits published since 2010, and GWAS was used to further validate the MQTL. Comparative genomics was used to identify wheat orthologs from rice and Arabidopsis thaliana to discover genomic regions and important candidate genes affecting flag leaf morphology in wheat. The aim of this study was to better understand the genetic mechanism underlying flag leaf morphological traits, and to provide useful information for genetic improvement of plant architecture and yield potential in wheat.

Integration of QTL for flag leaf morphological traits in wheat
In this study, the high-density map developed by Venske et al. [64] was used as the consensus map, which mainly includes three types of markers: SNP, DArT and SSR markers. SNP markers were derived from SNP array and genotyping-by-sequencing (GBS) [65,66]. SSR markers came from three genetic maps (Wheat, Consensus SSR 2004, Wheat Composite 2004 and

PLOS ONE
Wheat Synthetic × OPATA) in https://wheat.pw.usda.gov/GG3/ [67,68]. The diversity Array technology (DArT) marker was derived from the wheat consensus map 4.0 integrated by more than 100 genetic maps. According to the LOD value, phenotypic variation explained (PVE), confidence interval and position of QTL, the QTL for the target trait was mapped to the consensus map by using BioMercator v4.2 software [69], and the principle that the flanking marker of QTL interval corresponds to the consensus map interval was followed. Before mapping to the consensus map, the 95% confidence intervals (CI) of QTL identified in different studies were inferred by using the following formulas: (1) C.I. = 530 / (N×PVE); (2) C.I. = 163 / (N×PVE); (3) C.I. = 287 / (N×PVE), C.I. is the confidence interval of QTL, N is the size of mapping population, the value 530, 163 and 287 are specific population constants calculated by different simulations, formula (1), (2) and (3) is suitable for F 2 and backcross population, recombinant inbred line (RIL) population, and double haploid (DH) population, respectively [26,70]. Details of these initial QTLs are listed in S1 Table. QTL meta-analysis and verification by GWAS QTL meta-analysis for flag leaf morphological traits in wheat was carried out by using BioMercator v4.2 software. According to the number of QTL on each chromosome, two different analysis methods were used. When the QTL count on each chromosome is less than 10, MQTL is calculated for n independent QTLs by the method of Goffinet et al. [25]. Among the five models of 1, 2, 3, 4 and N, the lowest Akaike information criterion (AIC) value is considered as the best fitting model. When the QTL count on each chromosome exceeds 10, the method of Veyrieras et al. [71] is selected to determine the best QTL model based on AIC, AICc, AIC3, bayesian information criterion (BIC) and average weight of evidence (AWE), and the model with the lowest value of the selection criterion was used to determine MQTL. All the flanking markers sequences of MQTL were BLASTed against the wheat Chinese spring reference genome sequence (RefSeq v1.0) to obtain the physical position of MQTL. Ten papers published in the past five years on the genome-wide association studies of flag leaf morphology in wheat were collected (Table 2), and the physical location of the MTA (maker-traitassociation) in these studies was used to verify the accuracy of the MQTL region.

Mining candidate genes based on homology
According to the standard of mining highly reliable MQTLs by Venske et al. [71], MQTL with physical distance less than 20 Mb, genetic distance less than 1 cM and at least five overlapping QTLs were further selected as high confidence MQTL (hcMQTL). Combining with the genome annotation (https://wheat-urgi.versailles.inra.fr/seq-repository/annotations), the genes in the physical region of hcMQTL were analyzed, and the wheat orthologs in the physical region of hcMQTL were identified based on the genes related to flag leaf morphological traits of rice and Arabidopsis thaliana in Ensembl plant database (http://plants.ensembl.org/).

QTL integration for flag leaf morphological traits in wheat
A total of 465 QTLs related to flag leaf morphological traits were identified in the 20 papers published since 2010, covering 26 different mapping populations, among which 304 QTLs were projected into the consensus map. The number of QTL on each chromosome ranged from 2 on 3D to 38 on 5A, and the average number of QTLs on each chromosome was 14. Among them, 44.4% QTLs were distributed on A genome, 34.8% on B genome and 20.8% on D genome (Fig 1a). A total of 96 (31.6%), 104 (34.2%) and 80 QTLs (26.3%) were associated with flag leaf length, flag leaf width, and flag leaf area, respectively. Only 18 (5.9%) QTLs for flag leaf length-width ratio and 6 (2.0%) QTLs for flag leaf angle were identified (Fig 1b). The LOD score of individual QTLs ranged from 2.0 to 18.0, 54.28% of QTLs showed LOD scores from 3 to 4.5 (Fig 1c). The PVE range of individual QTL was 0.68-33.96%, and the PVE of 51.64% QTLs was within 3-9% (Fig 1d).

QTL meta-analysis for flag leaf morphological traits in wheat
A total of 304 QTLs were mapped to the consensus map, of which 275 QTLs were integrated into 55 MQTLs by meta-analysis, the remaining 29 QTLs were not integrated because they did not overlap with the above MQTLs (Table 3). These MQTLs were distributed on all chromosomes, with the number of MQTLs varying from one to four on each chromosome (Fig 2). The confidence interval of MQTL ranged from 0.06 to 16.45 cM, with the average interval size of 2.05 cM, which was 5.08-fold smaller than the initial QTL interval (Fig 3), the physical position interval ranged from 0.4 to 459.1 Mb, with the average physical distance of 66.5 Mb (Table 3).

MQTL validation by GWAS
Among the 55 MQTLs, 25 (45.45%) MQTLs were mapped into physical region smaller than 20 Mb in the wheat reference genome (Table 3). To determine the accuracy of MQTL region, GWAS studies on wheat flag leaf morphology reported in the past five years were used to verify MQTL. Since the decay distance of the wheat linkage disequilibrium was about 5 Mb, those overlapping MQTLs within 5 Mb of MTA were considered to be co-located with MQTL. Ten of the 25 MQTLs were verified in at least one GWAS study and co-located with 45 MTAs (Fig  4). The number of MTAs co-located in each MQTL ranged from 1 to 22, in which MQTL-10 co-located with 22 MTAs, followed by MQTL-13 that co-located with 11 MTAs. In addition, clusters or nested distributions of MQTL were observed, such as MQTL-6 (2A: 89.6-101.0

PLOS ONE
(S2 Table). In order to identify candidate genes related to leaf morphology among the four hcMQTLs, based on the Ensembl plant database (http://plants.ensembl.org/), 11 genes (six for rice and five for Arabidopsis) that regulate leaf morphology from the hcMQTL region were identified, of which, seven wheat orthologs for Osmtd1, three orthologs for FRS7, two orthologs for Roc8, and only one ortholog each for the other genes were found (Table 4).

Discussion
QTL meta-analysis developed by Goffinet et al.
[25] is a method for identifying consistent and stable QTLs and improving the accuracy of their genetic positions. The length, width, area, and angle of flag leaves are all important factors in determining wheat plant architecture and yield potential [5,[90][91][92]. Many genetic studies have been conducted to identify QTL for flag leaf morphological traits in wheat ( Table 1). Most of the initial QTLs collected in this study were distributed on A genome, and the least on D genome, which was slightly different from the results of previous studies regarding the distribution of initial QTLs for wheat yield and related traits on the genome (the most QTLs were distributed on B genome, but the least on D genome) [36,93], which might be due to the limited number of QTLs for flag leaf morphological traits, resulting in inconsistent results with previous studies. The less QTL on D genome may be related to the low-level polymorphism on D genome [94]. In this study, the maximum likelihood estimation method was used in meta-analysis in combination with the genetic locations of hundreds of QTLs for flag leaf morphological traits in wheat, and with consideration of population size and other QTL information, 275 of the 304 QTLs were mapped onto the consensus map and integrated into 55 MQTLs in wheat. Due to the pleiotropic effect of genes on flag leaf morphology in wheat, more than 90% (50/55) of the MQTLs were associated with at least two flag leaf morphological traits, and about 43.64% (24/55) of the MQTLs affected three or more flag leaf morphological traits simultaneously (Table 3).
After integrating QTLs by meta-analysis, the average confidence interval of MQTL was 2.05 cM, which was about 5.08-fold smaller than the average confidence interval (10.41 cM) of the initial individual QTL (Fig 3). Accordingly, the physical intervals of MQTL on chromosomes were further reduced, improving the accuracy of QTL mapping. The primary QTL mapping to fine mapping usually needs to increase molecular marker density [95] or construct fine mapping populations such as near-isogenic lines [96,97]. In certain cases, QTL meta-analysis could replace or enhance these approaches. For example, MQTL-52 was integrated by eight QTLs for flag leaf length, flag leaf width, and flag leaf area from two different populations and finally located within the interval of 58.64-59.52 cM on chromosome 7B, with the physical interval of 583.5-588.9 Mb, which was much smaller than the confidence interval of the initial QTL.
Compared with QTL linkage analysis mapping, linkage disequilibrium-based genome-wide association studies (GWAS) is another method for precisely locating genomic regions of quantitative traits. In previous studies, the results of wheat MQTL verification by GWAS have been reported [98,99]. For example, Aduragbemi et al. [100] identified 51 MTA and 29 MQTLs colocated for leaf rust resistance loci using GWAS. Yang et al. [101] verified MQTL for wheat yield and yield-related traits using GWAS results published in recent years, and found that about 60% of MQTLs were co-located with MTA. In this study, based on the GWAS results of wheat flag leaf morphological traits published in recent years, 45 MTAs and 10 MQTLs were identified, which indicated that these genomic regions controlling flag leaf morphological traits might be less affected by the genetic background and environment. The 10 MQTLs verified by GWAS provided a basis for the accurate mining candidate genes that affect flag leaf morphology in wheat. Loffler et al. [102] proposed the criteria for selection of MQTL for use in breeding programs: the MQTL with confidence interval genetic distance less than 2 cM, no less than 4 initial QTLs from different studies with PVE> 10%. On this basis, we determined three potential MQTLs, MQTL-1, MQTL-13 and MQTL-25, that could be used to improve wheat flag leaf morphological traits.

PLOS ONE
Flag leaf morphology is one of the important traits of plant architecture in wheat breeding. Moreover, previous studies reported the correlation between flag leaf morphology and plant structure traits such as plant height and tiller number [57,78]. Hu et al. [57] revealed the genetic mechanism of yield-related traits in wheat using four RIL populations and found that plant height had a significant positive correlation with FLL and FLW, and QTL affecting both plant height and FLW were detected at 0-3.5 cM on chromosome 5A. In addition, Muhammad et al. [78] identified five SNP markers affecting PH, FLL and FLW simultaneously on chromosomes 1A, 3A, 3B, 5A, and 6B in natural populations of wheat. Some previously reported major QTLs and genes controlling plant height and tillering number in wheat were identified in the hcMQTL region in this study. The gene Csl-1A (chr1A:6.4 Mb) controlling the tiller number in wheat [103] was identified near MQTL-1. Saini et al. These examples indicate that these three hcMQTLs may carry some major genes that improve the plant architecture of wheat, such as plant height and tiller.

PLOS ONE
In cereal with large complex genomes, such as wheat, barley and maize, localization based on homologous cloning is an effective way to identify important genes associated with complex traits. With the wide application of high-throughput sequencing technology, many crops genome sequences have been published, which conduces to identify conserved genome regions and key genes in different crops. For example, the rice OsLG1 gene encodes a SBP DNA binding protein, which affects the development of auricle and ligule [104], and the orthologous gene TaSPL8 in wheat has also been found to have similar function in rice [23].

PLOS ONE
In this study, a total of 20 wheat orthologs were identified in four hcMQTLs, 10 of which were low-confidence, and the annotation information might be inaccurate. The remaining 10 wheat orthologs were potential candidate genes for regulating the leaf morphology in wheat, including 7 from rice genes and 3 from Arabidopsis genes (Table 4). In total, 7 wheat orthologs of the rice gene Osmtd1 were located in the MQTL-1 region, namely TraesCS1A02G022100, TraesCS1A02G022200, TraesCS1A02G022300, TraesCS1A02G024800, TraesCS1A02G024900, TraesCS1A02G025000 and TraesCS1A02G025300. Osmtd1 gene encodes a putative inhibitor I family protein regulating rice tillering and leaf angle [80]. Therefore, these seven wheat orthologs may be reliable candidates for regulating wheat leaf angle as the Osmtd1 gene in rice.
The wheat ortholog of Arabidopsis PLL5 gene, TraesCS1A02G033700, is located in the MQTL-1 region and encodes a protein belonging to the phosphatase 2C family, which regulates leaf development. The mutant pll5 has shorter, narrower and curlier leaves than the wildtype leaves [84]. Hence, it suggested that TraesCS1A02G033700 is a credible candidate gene affecting leaf development in wheat. The CAND1 gene encodes unmodified CUL1-interacting protein in Arabidopsis, and participates in many developmental pathways controlled by ubiquitin/proteasome-mediated degradation of protein [87]. The rosette leaves of cand1 mutants are much smaller than that of wild-type plants and have a wavy morphology. The ortholog TraesCS2B02G051000 of wheat located in MQTL-11 region encodes CUL1-related NEDD8 dissociation protein, which might be a candidate gene affecting the wheat leaves morphology. The ECT3 gene encodes the YTH domain protein in Arabidopsis, which has been previously proved to be related to leaf morphogenesis in Arabidopsis [89]. The wheat ortholog TraesCS2D02G012200, located in the MQTL-13 region, might be a reliable candidate gene involved in the regulation of leaf development.
In conclusion, using the high-density integration map developed by Venske et al. [64] as the consensus map and QTL meta-analysis, we integrated the QTL for flag leaf morphological traits previously identified in wheat, and validated 10 MQTLs with GWAS information. Three potential MQTLs, MQTL-1, MQTL-13 and MQTL-25 that regulate flag leaf morphological traits were identified in this study. These MQTL flanking markers can be used for molecular marker assisted breeding to improve flag leaf morphological traits in wheat. Furthermore, using functional annotation information from genes within the hcMQTL interval and a comparative genomics strategy, ten wheat orthologs were identified as potential candidate genes affecting wheat flag leaf morphology, providing potential targets for fine mapping, and gene cloning.
Supporting information S1